#!/usr/bin/python
# plot_SimpleIntegrate.py

#import math
import numpy
from matplotlib import use as useBackend
useBackend('Agg')
import matplotlib.pyplot as ppt


def cheat111():
    '''Cheat to clear the plot.'''
    ppt.subplot(311)
    ppt.plot([0,1],[1,2])
    ppt.subplot(131)
    ppt.plot([0,1],[1,2])
    ppt.subplot(111)
    ppt.hold(True)


colorlist = ['k', 'b', 'g', 'r', 'c', 'm', 'y', 'indigo', 'darkorange', 'saddlebrown', 'forestgreen', 'deeppink', 'slategray', 'maroon', 'lime', 'bisque', 'aquamarine', 'darkcyan', 'darkmagenta', 'mediumpurple', 'olive', 'purple', 'salmon'] #23


def choosecolor(j, colorlist):
    ncolors = len(colorlist)
    while True:
        if j > (ncolors - 1):
            j = j - ncolors
        else:
            colornow = colorlist[j]
            break
    return colornow

def safeindex(x,y, options=None):
    whacked = False
    if len(x) != len(y):
        print "You are wrong.", x.shape, y.shape
        print "I am deliberately crashing this because you screwed up."
        whacked = True
    imin = 0
    imax = len(x)
    for i in range(0, len(x)):
        if x[i] == 0.0:
            imin = int(i)
        else:
            break
    for i in range(0, len(x)):
        thisi = len(x) - i - 1
        if x[thisi] == 0.0:
            imax = int(thisi)
        else:
            break
    if imin == (len(x) - 1):
        imin = imin + 2
        imax = imin + 1
    jmin = 0
    jmax = len(y)
    for j in range(0, len(y)):
        if y[j] == 0.0:
            jmin = int(j)
        else:
            break
    for j in range(0, len(y)):
        thisj = len(y) - j - 1
        if y[thisj] == 0.0:
            jmax = int(thisj)
        else:
            break
    if jmin == (len(y) - 1):
        jmin = jmin + 2
        jmax = jmin + 1
    if options==None:
        trumin = max(imin,jmin)
        trumax = min(imax,jmax)
    if options=='y':
        trumin = jmin
        trumax = jmax
    if options=='x':
        trumin = imin
        trumax = imax
    if whacked == True:
        trumax = 'wrench'
    return (trumin, trumax)
    

HanMvLims = [6, 18.25, 3e-4, 6]
HanDSLims = [0, 2, 3e-4, 1]
HanDLLims = [0, 2, 3e-4, 1]
HanMassLims = [7e-2, 2., 3e-4, 1]

def plot(simpars, phypars, ratestuple, UseHanLims=False):
    '''Shows GM_SimpleIntegrate's results.'''
    for i in simpars:
        cmd = "%s = simpars['%s']" % (i,i)
        exec cmd
    for i in phypars:
        cmd = "%s = phypars['%s']" % (i,i)
        exec cmd
    print "plot_SimpleIntegrate is going about its business."
    (lensrates, MvHist, MassHist, DSHist, DLHist, bHist) = ratestuple
    # Set up naming/identifying string and plot legends
    sourcestring = ''
    legendlist = []
    plotlist = []
    for i in range(0, len(SourceMag)):
        sourcestring = sourcestring + '_' + str(SourceMag[i])
        legendname = 'Mv = ' + str(SourceMag[i])
        legendlist.append(legendname)
        cn = choosecolor(i, colorlist)
        thisplot = ppt.plot([1,2], [0,1], c=cn, label=legendname)
        plotlist.append(thisplot)
    if UseHanLims == True:
        sourcestring = '_HanLims' + sourcestring
    # Rate contribution as a function of mass
    cheat111()
    masses = MvMass[:,1]
    for i in range(0, len(MassHist)):
        cn = choosecolor(i, colorlist)
        (trumin, trumax) = safeindex(masses, MassHist[i])
        ppt.loglog(masses[trumin:trumax], MassHist[i][trumin:trumax], c=cn)
        ppt.loglog(masses[trumin:trumax], MassHist[i][trumin:trumax], c='k', ls='none', marker=',')
    if UseHanLims == True:
        ppt.axis(HanMassLims)
#    ppt.legend(plotlist, legendlist)
    ppt.legend()
    ppt.xlabel('Mass in Msun')
    ppt.ylabel('Contribution to lensing rate per mass interval')
    ppt.title('Lensing rate contribution as a function of mass')
    filename = outdir+'rates_masses'+sourcestring+'.pdf'
    ppt.savefig(filename)
    # Rate contribution as a function of source apparent magnitude
    cheat111()
    for i in range(0, len(MvHist)):
        lenMv = len(MvHist[i])
#        Mvxvals = limitset[:lenMv,:].mean(axis=1)
        Mvxvals = numpy.array(MvVals)
        cn = choosecolor(i, colorlist)
        (trumin, trumax) = safeindex(Mvxvals, MvHist[i], 'y')
        ppt.semilogy(Mvxvals[trumin:trumax], MvHist[i][trumin:trumax], c=cn)
        ppt.semilogy(Mvxvals[trumin:trumax], MvHist[i][trumin:trumax], c='k', ls='none', marker=',')
    if UseHanLims == True:
        ppt.axis(HanMvLims)
    ppt.legend(plotlist, legendlist)
    ppt.xlabel('Unlensed apparent source magnitude')
    ppt.ylabel('Contribution to lensing rate per Mv interval')
    ppt.title('Lensing rate contribution as a function of source Mv')
    filename = outdir+'rates_Mv'+sourcestring+'.pdf'
    ppt.savefig(filename)
    # Rate contribution as a function of DS
    cheat111()
    Distance = dlocs*PtScaled/1000. # pc -> kpc
    ###DSHist = DSHist*1000 # pc -> kpc # Unnecessary b/c normalized?
    for i in range(0, len(DSHist)):
        cn = choosecolor(i, colorlist)
        (trumin, trumax) = safeindex(Distance, DSHist[i], 'y')
        ppt.semilogy(Distance[trumin:trumax], DSHist[i][trumin:trumax], c=cn)
        ppt.semilogy(Distance[trumin:trumax], DSHist[i][trumin:trumax], c='k', ls='none', marker=',')
    if UseHanLims == True:
        ppt.axis(HanDSLims)
    ppt.legend(plotlist, legendlist)
    ppt.xlabel('Distance to source, kpc')
    ppt.ylabel('Contribution to lensing rate per distance interval')
    ppt.title('Lensing rate contribution as a function of source distance')
    filename = outdir+'rates_DS'+sourcestring+'.pdf'
    ppt.savefig(filename)
    # Rate contribution as a function of DL
    cheat111()
#    Distance = dlocs*PtScaled/1000. # pc -> kpc
    ###DLHist = DLHist*1000 # pc -> kpc # Unnecessary b/c normalized?
    for i in range(0, len(DLHist)):
        cn = choosecolor(i, colorlist)
        (trumin, trumax) = safeindex(Distance, DLHist[i], 'y')
        ppt.semilogy(Distance[trumin:trumax], DLHist[i][trumin:trumax], c=cn)
        ppt.semilogy(Distance[trumin:trumax], DLHist[i][trumin:trumax], c='k', ls='none', marker=',')
    if UseHanLims == True:
        ppt.axis(HanDLLims)
    ppt.legend(plotlist, legendlist)
    ppt.xlabel('Distance to lens, kpc')
    ppt.ylabel('Contribution to lensing rate per distance interval')
    ppt.title('Lensing rate contribution as a function of lens distance')
    filename = outdir+'rates_DL'+sourcestring+'.pdf'
    ppt.savefig(filename)
    print "plot_SimpleIntegrate has concluded its business."
